Reduced-form VAR Modeling
==============================


Description
-------------

.. math::
	
	y_{t}=C\left( r_{t}\right) x_{t}+B_{1}\left( r_{t}\right)
	y_{t-1}+...+B_{p}\left( r_{t}\right) y_{t-p}+u_{t}

With :math:`r_{t}=1,2,...,h` and transition probabilities :math:`p_{r_{t},r_{t+1}}\left( I_{t}\right)`

Quick-start examples
----------------------

A constant-parameter VAR
~~~~~~~~~~~~~~~~~~~~~~~~~~

Collecting and transforming data
"""""""""""""""""""""""""""""""""""
::

	% CLVMNACSCAB1GQNO : GDP Norway
	% IR3TIB01NOQ156N : 3-month Interbank interest rate
	% NORCPGRLE01IXOBQ : CPI excluding food and energy
	% CCUSSP01NOQ650N : Spot Exchange rate : 1NOK = x USD
	% POILWTIUSDQ : Global price of WTI Crude
	xrange='1990Q1:2022Q3';
	rawdb=fetch_fred({'NORCPGRLE01IXOBQ','IR3TIB01NOQ156N','CLVMNACSCAB1GQNO',...
	    'CCUSSP01NOQ650N','POILWTIUSDQ'});
	rawdb1.P=rawdb(1).series(xrange);
	rawdb1.INTRATE=rawdb(2).series(xrange);
	rawdb1.Y=rawdb(3).series(xrange);
	rawdb1.EXRATE=rawdb(4).series(xrange);
	rawdb1.POIL=rawdb(5).series(xrange);
	rawdb=rawdb1; clear rawdb1
	
	db=struct();
	db.PAI=rawdb.P/lag(rawdb.P,1);
	db.R=1+rawdb.INTRATE/100;
	db.GROWTH=rawdb.Y/lag(rawdb.Y,1);
	db.EXRATE=1/rawdb.EXRATE; % so that increase = depreciation of NOK
	db.PAIOIL=rawdb.POIL/lag(rawdb.POIL,1);

Setting up the Reduced.form VAR
""""""""""""""""""""""""""""""""""
::

	endog={'PAIOIL','GROWTH','PAI','R','EXRATE'};
	
	exog={};
	
	nlags=4;
	
	const=true;
	
	mdl = rfvar(endog, exog, nlags, const);

Estimating the VAR using classical techniques
"""""""""""""""""""""""""""""""""""""""""""""""""" 
::

	data_range={db.GROWTH.start,db.GROWTH.finish};
	
	mdlest=estimate(mdl,db,data_range);

Restrictions on the VAR : Domestic variables do not affect oil prices
""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""

::

	linres={
	    'b1(PAIOIL,PAI)=0'
	    'b1(PAIOIL,GROWTH)=0'
	    'b1(PAIOIL,R)=0'
	    'b1(PAIOIL,EXRATE)=0'
	    'b2(PAIOIL,PAI)=0'
	    'b2(PAIOIL,GROWTH)=0'
	    'b2(PAIOIL,R)=0'
	    'b2(PAIOIL,EXRATE)=0'
	    };
	
	% or programmatically
	% linres=cell(0,1);
	% for ilag=1:nlags
	%     for iv=2:numel(endog)
	%         y=endog{iv};
	%         linres{end+1,1}=sprintf('b%0.0f(PAIOIL,%s)=0',ilag,y); %#ok<SAGROW>
	%     end
	% end

Estimate the restricted VAR
"""""""""""""""""""""""""""""""
::

	prior=[];
	
	mdlest_restr=estimate(mdl,db,data_range,prior,linres);

Identification
"""""""""""""""""
::

	shock_names={'oilp','demand','costpush','mp','forex'};
	
	ident_restr1={
	    % normalization with sign restrictions
	    'PAIOIL{0}@oilp','+'
	    'GROWTH{0}@demand','+'
	    'PAI{0}@costpush','+'
	    'R{0}@mp','+'
	    'EXRATE{0}@forex','+'
	    % first set
	    'PAIOIL{0}@demand',0
	    'PAIOIL{0}@costpush',0
	    'PAIOIL{0}@mp',0
	    'PAIOIL{0}@forex',0
	    % second set
	    'GROWTH{0}@costpush',0
	    'GROWTH{0}@mp',0
	    'GROWTH{0}@forex',0
	    % third set
	    'PAI{0}@mp',0
	    'PAI{0}@forex',0
	    % fourth set
	    'R{0}@forex',0
	    };
	
	agnostic=true;
	
	max_trials=6000;
	
	Rfunc=struct(); ident=struct();
	
	[Rfunc.unrestr,ident.unrestr]=identification(mdlest,ident_restr1,shock_names,...
	    agnostic,max_trials);
	
	[Rfunc.restr,ident.restr]=identification(mdlest_restr,ident_restr1,shock_names,...
	    agnostic,max_trials);
	
	disp(ident.restr)
	
	disp(ident.unrestr)

Compare structural shocks
"""""""""""""""""""""""""""""
::

	params=[];
	sshocks=struct();
	sshocks.unrestr=structural_shocks(mdlest,...
	    params,Rfunc.unrestr,shock_names);
	sshocks.restr=structural_shocks(mdlest_restr,...
	    params,Rfunc.restr,shock_names);
	
	% plots
	%-------
	
	titel='Structural shocks';
	
	figure('name',titel);
	
	for ii=1:numel(shock_names)
	    
	    thisname=shock_names{ii};
	    
	    subplot(3,2,ii)
	    
	    d=[sshocks.unrestr.(thisname),sshocks.restr.(thisname)];
	    
	    plot(d,'linewidth',2)
	    
	    xtickangle(45)
	    
	    title(thisname)
	    
	    if ii==1
	        
	        legend({'unrestr','restr'})
	        
	    end
	    
	end
	
	% xrotate(45)
	
	[~,h]=sup_label(titel,'t');
	
	set(h,'fontsize',12)

Compare CHOLESKI irfs
"""""""""""""""""""""""""
::

	cholShocks=[];
	
	myirfs=irf([mdlest,mdlest_restr],cholShocks,40);
	
	cholShocks=fieldnames(myirfs);
	
	tex=struct();
	tex.PAIOIL='Oil price inflation';
	tex.GROWTH='GDP growth';
	tex.PAI='Inflation';
	tex.R = 'Interest rate';
	tex.EXRATE='exchange rate';
	
	for ishock=1:numel(cholShocks)
	    
	    shock=cholShocks{ishock};
	    
	    titel=['Impulse responses to a ',shock,' shock'];
	    
	    figure('name',titel);
	    
	    for ii=1:numel(endog)
	        
	        subplot(3,2,ii)
	        
	        plot(myirfs.(shock).(endog{ii}),'linewidth',2)
	        
	        title(tex.(endog{ii}))
	        
	        if ii==1
	            
	            legend({'unrestr','restr'})
	            
	        end
	        
	    end
	    
	    [~,h]=sup_label(titel,'t');
	    
	    set(h,'fontsize',12)
	    
	end	

Compare irfs based on the identification scheme
"""""""""""""""""""""""""""""""""""""""""""""""""""
::

	params=[];
	
	% Note the two models have the same identification !!!
	% This is why we can run them together
	myirfs=irf([mdlest,mdlest_restr],shock_names,40,params,Rfunc.unrestr);
	
	for ishock=1:numel(shock_names)
	    
	    shock=shock_names{ishock};
	    
	    titel=['Impulse responses to a ',shock,' shock'];
	    
	    figure('name',titel);
	    
	    for ii=1:numel(endog)
	        
	        subplot(3,2,ii)
	        
	        plot(myirfs.(shock).(endog{ii}),'linewidth',2)
	        
	        title(tex.(endog{ii}))
	        
	        if ii==1
	            
	            legend({'unrestr','restr'})
	            
	        end
	        
	    end
	    
	    [~,h]=sup_label(titel,'t');
	    
	    set(h,'fontsize',12)
	    
	end

Variance decomposition
"""""""""""""""""""""""
::

	params=[];
	
	vd=variance_decomposition(mdlest_restr,params,Rfunc.restr);
	
	% plot decompositions
	%---------------------
	
	range='0:50'; % pick a range for the plots
	
	figure('name','Variance Decomposition');
	
	for iv=1:numel(endog)
	    
	    d=vd.conditional.(endog{iv});
	    
	    subplot(3,2,iv)
	    
	    plot_decomp(range,d)
	    
	    if iv==1
	        
	        legend(shock_names,'location','SE',...
	            'Orientation','horizontal')
	        
	    end
	    
	    title(tex.(endog{iv}))
	    
	end

Historical decomposition
""""""""""""""""""""""""""""
::

	params=[];
	
	hd=historical_decomposition(mdlest_restr,params,Rfunc.restr);
	
	% plot decompositions
	%----------------------
	
	shock_only=true;
	
	titel='Model with block exogeneity: Historical Decomposition';
	
	figure('name',titel);
	
	for iv=1:numel(endog)
	    
	    d=hd.(endog{iv});
	    
	    subplot(3,2,iv)
	    
	    if shock_only
	        
	        shock_tex=shock_names;
	        
	        d=d(shock_names);
	        
	    else
	        
	        shock_tex=d.varnames;
	        
	    end
	    
	    plot_decomp(d)
	    
	    if iv==1
	        
	        legend(shock_tex,'location','SE','Orientation','horizontal')
	        
	    end
	    
	    title(tex.(endog{iv}))
	    
	end
	
	[~,h]=sup_label(titel,'t');
	
	set(h,'fontsize',12)

Bootstrap
""""""""""""""
::

	n=1000;
	    
	params=bootstrap(mdlest_restr,n);

Variance decomposition distribution
"""""""""""""""""""""""""""""""""""""""
::

	ci=[30,50,68,90];
	
	vd=variance_decomposition(mdlest_restr,params,Rfunc.restr);
	
	shock_tex=shock_names;
	
	myrange='1:50';
	
	for iv=1:numel(endog)
	    
	    vname=tex.(endog{iv});
	    
	    titel=['Variance Decomposition (in %) of ',vname];
	    
	    figure('name',titel);
	    
	    d=vd.conditional.(endog{iv});
	    
	    d.varnames=shock_names;
	    
	    d=pages2struct(d);
	    
	    contributors=fieldnames(d); % = shock_names
	    
	    for ii=1:numel(contributors)
	    
	        subplot(3,2,ii)
	        
	        % note we are multiplying by 100, this just by pure convenience
	        %--------------------------------------------------------------
	        out=fanchart(100*d.(contributors{ii})(myrange),ci);
	        
	        plot_fanchart(out)
	        
	        title(contributors{ii})
	        
	        axis tight
	        
	    end
	    
	    [~,h]=sup_label(titel,'t');
	    
	    set(h,'fontsize',12)
	    
	end

Historical decomposition distribution
""""""""""""""""""""""""""""""""""""""""""
::

	ci=[30,50,68,90];
	
	% compute decompositions
	    
	    hd=historical_decomposition(mdlest_restr,...
	        params,Rfunc.restr);
	
	% plot decompositions
	
	for iv=1:numel(endog)
	    
	    vname=tex.(endog{iv});
	    
	    titel=['Model with restrictions : Historical Decomposition of ',vname];
	    
	    figure('name',titel);
	    
	    d=pages2struct(hd.(endog{iv}));
	    
	    contributors=fieldnames(d);
	    
	    for ii=1:numel(contributors)
	        
	        subplot(5,2,ii)
	        
	        out=fanchart(d.(contributors{ii}),ci);
	        
	        plot_fanchart(out)
	        
	        xtickangle(45)
	        
	        title(contributors{ii})
	        
	        axis tight
	        
	    end
	    
	    [~,h]=sup_label(titel,'t');
	    
	    set(h,'fontsize',12)
	    
	end

IRF distribution
"""""""""""""""""""""
::

	    myirfs=irf(mdlest_restr,shock_names,...
	        40,params,Rfunc.restr);
	    
	% IRFs plots
	%--------------
	    
	for ishock=1:numel(shock_names)
	    
	    shock=shock_names{ishock};
	    
	    titel=['Model with restrictions : IRFs to a ',shock,' shock'];
	    
	    figure('name',titel);
	    
	    for ii=1:numel(endog)
	        
	        subplot(3,2,ii)
	        
	        d=myirfs.(shock).(endog{ii});
	        
	        out=fanchart(d,ci);
	        
	        plot_fanchart(out)
	        
	        title(tex.(endog{ii}))
	        
	        axis tight
	        
	    end
	    
	    [~,h]=sup_label(titel,'t');
	    
	    set(h,'fontsize',12)
	    
	end

Bayesian estimation
"""""""""""""""""""""""
::

	v=rfvar(endog,exog,nlags,const);
	
	%set prior
	
	var_prior=rfvar.prior_template();
	% modify as needed
	
	var_prior.type='sz';
	
	prior=struct('var',var_prior);
	% prior.nonvar
	% prior.endogenous
	
	% unrestricted model
	%--------------------
	
	ve=estimate(v,db,data_range,prior);%,restrictions
	
	% restricted model
	%------------------
	
	ve_lr=estimate(v,db,data_range,prior,linres);

Posterior sampling of parameters
""""""""""""""""""""""""""""""""""""""
::

	params=struct();
	params.ve=ve.estim_.sampler(1000);
	params.ve_lr=ve_lr.estim_.sampler(1000);

Bayesian forecasting
""""""""""""""""""""""
::

	myfkst=struct();
	date_start='2003Q1';
	myfkst.ve=forecast(ve,db,date_start,params.ve);
	
	myfkst.ve_lr=forecast(ve_lr,db,date_start,params.ve_lr);
	
	% set environment
	
	ci=[30,50,68,90];
	
	modelnames={'ve','ve_lr'};
	
	% Forecast plots
	
	for jj=1:numel(modelnames)
	    
	    modname=modelnames{jj};
	    
	    figure('name',['model (',modname,') Forecasts of Norwegian Data']);
	    
	    for ii=1:numel(endog)
	        
	        subplot(3,2,ii)
	        
	        d=myfkst.(modname).(endog{ii});
	        
	        out=fanchart(d,ci);
	        
	        plot_fanchart(out,[244, 122, 66]/244)
	        
	        title(tex.(endog{ii}))
	        
	        axis tight
	        
	    end
	    
	    xrotate(45)
	    
	end

Conditional forecasting
""""""""""""""""""""""""""""
::

	myfkst=struct();
	% date_start=[];
	date_start='2003Q1';
	nsteps=12;
	shock_uncertainty=false;
	Rfunc=[]; % No need for identification
	
	conditions=struct();
	conditions.R={'2003Q1','2004Q4'}; % range over which we want to condition
	
	myfkst.ve=forecast(ve,db,date_start,params.ve,nsteps,...
	    shock_uncertainty,Rfunc,conditions);
	
	myfkst.ve_lr=forecast(ve_lr,db,date_start,params.ve_lr,nsteps,...
	    shock_uncertainty,Rfunc,conditions);
	
	% set environment
	
	ci=[30,50,68,90];
	
	% Forecast plots
	
	for jj=1:numel(modelnames)
	    
	    modname=modelnames{jj};
	    
	    figure('name',['model (',modname,') Forecasts of Norwegian Data']);
	    
	    for ii=1:numel(endog)
	        
	        subplot(3,2,ii)
	        
	        d=myfkst.(modname).(endog{ii});
	        
	        out=fanchart(d,ci);
	        
	        plot_fanchart(out,[244, 122, 66]/255)
	        
	        title(tex.(endog{ii}))
	        
	        axis tight
	        
	    end
	    
	    xrotate(45)
	    
	end

Adding regime switching with constant transition probabilities
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
::

	%% adding constant transition probabilities
	
	mc=struct();
	mc.name='policy';
	mc.number_of_states=2;
	mc.controlled_parameters={'b(R,:)'};
	mc.endogenous_probabilities=[];
	mc.probability_parameters=[];
	
	mdl = rfvar(endog, exog, nlags, const,mc);

::

	%% Set priors and Estimate the VAR
	
	switch_prior=struct();
	switch_prior.policy_tp_1_2={0.5,0.1,0.3,'beta'};
	switch_prior.policy_tp_2_1={0.5,0.1,0.3,'beta'};
	
	prior=struct();
	prior.nonvar=switch_prior;
	
	% ___Optional prior on the VAR__
	% var_prior=svar.prior_template();
	% var_prior.type='sz';
	% prior.var=var_prior;
	
	data_range={db.GROWTH.start,db.GROWTH.finish};
	
	mdlest=estimate(mdl,db,data_range,prior);


Adding regime switching with time~varying transition probabilities
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
::

	%% adding time-varying transition probabilities
	
	mc=struct();
	mc.name='policy';
	mc.number_of_states=2;
	mc.controlled_parameters={'b(R,:)'};
	mc.endogenous_probabilities={
	    'policy_tp_1_2=1/(1+exp(a12*(PAI-c12)))'
	    'policy_tp_2_1=1/(1+exp(a21*(PAI-c21)))'
	    };
	mc.probability_parameters={'a12','a21','c12','c21'};
	
	mdl = rfvar(endog, exog, nlags, const,mc);

::

	%% Set priors and estimate the RFVAR
	
	switch_prior=struct();
	switch_prior.a12={0,0,1,'normal'};
	switch_prior.a21={0,0,1,'normal'};
	switch_prior.c12={1.05^.25,1.05^.25,0.05,'gamma'};
	switch_prior.c21={1.05^.25,1.05^.25,0.05,'gamma'};
	
	prior=struct();
	prior.nonvar=switch_prior;
	
	% ___Optional prior on the VAR__
	% var_prior=svar.prior_template();
	% var_prior.type='sz';
	% prior.var=var_prior;
	
	data_range={db.GROWTH.start,db.GROWTH.finish};
	
	mdlest=estimate(mdl,db,data_range,prior);

Technical documentation for rfvar objects
-------------------------------------------

.. toctree::
   :maxdepth: 2
   :caption: Contents:

   rfvar_properties_methods